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I. INTRODUCTION 


During the last few years impressive progress was made 
in the experimental in vest igation of ultracold atomic 
gases in optical lattices.^ It is now possible to mea¬ 
sure not only density profiles and static correlation func¬ 
tions but even dynamical quantities such as spectral func¬ 
tions of the trapped particles, using Bragg spectroscopy.l^ 
Thereby momentum resolved bosonic spectral functions 
were obtained in the case of condensed boson^^H^ and 
across the phase transition to a Mott insulator.l^ By 
contrast, only rather few calculations of bosonic spec¬ 
tral functions have been performed so far for the Bose- 
Hubbard model. Early studies were based on the strong 
coupling approximation to the Bose-Hubbard model.lSl 
Later the weak-coupling limit at zero temperature was 
thoroughly analyzed with the functional renormaliza¬ 
tion group Other approaches include the variational 
cluster approximation, 1131 which was used to investigate 
systems in on^i^ and twcP^Hfil dimensions, the quan¬ 
tum rotor approximation,which was recently applied 
to three-dimensional systems ,13^ and the linked-cluster 
expansion (LCE).li^ In particular the LCE, which was 
employed to approximately solve the equations of the 
bosonic dynamical mean-field theory (B-DMET),1^ al¬ 
lows one to study the Bose-Hubbard model for strong 
interactions and near the phase boundary between the 
superfluid and the Mott-insulating phase. The B-DMFT 
treats local correlations in time exactly and includes spa¬ 
tial correlations on the mean-field level. It is applicable 
for all values of the Hubbard interaction, density and 
temperature. The B-DMFT was de rived subsequent to 
the DMFT for lattice fermMn^^illll and was appli ed to 
various bosonic problem^^lMl as well as to mixture^^HHl 
of bosons and fermions. It has been extended and ap¬ 
plied also to inhomogeneous situation^^ and to bosonic 
systems in non-equilibriunP^. 


In this paper we solve the Bose-Hubbard model within 
the B-DMFT framework on a simple cubic lattice in three 
dimensions using a continuous-time quantum Monte 
Carlo (CT-QMC) solver.!^ This approach is known to 
give exce llent results for the phase diagram and static 
propertie^22Iini of interacting bosons. We compute the 
momentum resolved and the momentum integrated spec¬ 
tral functions and the dispersion relation of interacting 
bosons in, both, the superfluid and the Mott insulating 
phase. The following questions will be addressed and 
answered: (i) How well can strong-coupling approaches, 
which are known to capture the phase diagram of cor¬ 
related bosons very well, describe dynamical properties 
such as spectral functions? (ii) How does the presence of 
the superfluid influence the spectral properties of normal 
bosons? (iii) How is the dispersion relation modified by 
the interaction? 


The paper is organized as follows: In Section |TT] we 
introduce the Bose-Hubbard model and recapitulate the 
main steps of the B-DMFT framework. We also intro¬ 
duce an improved method for calculating the self-energy 
which makes use of two-particle Green functions and dis¬ 
cuss the method for the numerical analytic continuation. 
In Section HI we present the results for the momentum 
resolved spectral function in the limit of weak, interme¬ 
diate, and strong interactions, respectively. Finally, Sec¬ 
tion ITV] concludes the paper with a summary. 


II. MODEL AND INVESTIGATION METHOD 


We consider spinless bosons on a lattice described by 
the Bose-Hubbard Hamiltonian 

^ = - E f E (1) 
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where b] (bi) is a bosonic creation (annihilation) oper¬ 
ator on a lattice site i, fj, is the chemical potential, U 
is the local interaction strength, and tij is the hopping 
amplitude. We assume nearest neighbour (NN) hopping, 
i.e., tij = t > 0 if sites i and j are NN and 0 other¬ 
wise. Our calculations were performed for a simple cubic 
lattice with coordination number z = 6 and t = 0.5 in 
arbitrary units. 


A. The B-DMFT action 


In the B-DMFT the lattice problem is replaced by 
a single-site (“impurity”) problem with self-consistency 
conditions. A detailed derivation can be found in Ref. 
[201 The impurity action reads 
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Here /3 = 1/T is the inverse of the temperature T, k = 
tij = zt is a geometrical parameter depending on a 
lattice type, and r is the imaginary (Matsubara) time. 
We use the Nambu vector notation for operators b and 
complex variables 5, i.e.. 



Then the connected Green functions are defined by 

G(t) = -(r,b(r)b^(0))^ 

= -(r.b(r)bt(0))-f (b(r))(bt(0)). 


ity at site 0 reads 

’4f=(b)(°). (6) 

The connection between the condensate field on a full 
lattice, $ = (b), and the condensate field on a lattice 
with a cavity, ’S', is expressed by 


’S' = 


drA(r) 




(7) 


Once the mean fields ^ and A(t) are known it is pos¬ 
sible to solve the impurity problem and calculate the lo¬ 
cal Green functions and local condensa te fie lds. This 
can be achieved by exact dia gonali zatiorP^*^, the LGE 
approximatiorfi^, or GT-QMG.I^^Ml jn this paper we use 
the GT-QMG method to solve the impurity problem. Our 
code is based on the work of Anders et We employ an 

improved calculation of the self-energy via two-particle 
Green functions, to be discussed in Section IIG and in 
the Appendix. Tests of our code show numerical agree¬ 
ment with the results of Ref. [231 


B. Self-consistency 


Once the single-particle Green functions have been ob¬ 
tained by solving the single-site problem, the local self¬ 
energy can be calculated from the Dyson equation 

I(ia.„) = ^ ^ - A(i.;J - [G(iccJ]-i, 

( 8 ) 

where LOn are even Matsubara frequencies for bosons 
{ujn = 27rn//3). In the B-DMFT the self-energies are 
local, i.e., momentum independent. The momentum re¬ 
solved Green functions are then given by 


where (...) denotes the equilibrium average in the grand 
canonical ensemble, and G(/3-|-r) = G(T).*^The last two 
terms in Eq. § represent the coupling of the site to two 
types of external mean-fields: (i) the Bose-Einstein con¬ 
densate (BEG) which is represented by a static mean-field 
and (ii) the dynamical mean-field of normal bosons 
represented by the matrix the elements of which 

are hybridization functions. 

The hybridization functions are related to the Green 
functions G*^°^ for a lattice with a cavity at site 0 (i.e., 
where site 0 is removed) bjESl 

A{t -t')= tiotjoGf^T - t'). (5) 

Similarly, the condensate field ’®' for a lattice with a cav- 


G(k,ia;„) = 


f iuJn + fJ, - ey: 

0 ^ 

0 
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where Ck is the dispersion relation for noninteracting 
bosons. To close the self-consistency equations we cal¬ 
culate the local Green functions from the k-integrated 
Dyson equation 


G(ia;„) = Y = 


iuJn + — € 0 

0 -ioJn + Ai - e 


— itnUJr, 


( 10 ) 

where D{e) is the noninteracting density of states (DOS). 
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particle, disconnected Green functions F(r) given by 

F(r) = -/r. f [6t6S](r)6(0) \ \ 

\ ^ y6t6tS](^)6t(o) lb^b^](T)b(0)J /' 

Details of this derivation are presented in Appendix 
In Fig. [^we see that the improved method significantly 
reduces the error, i.e., it now grows only linearly in fre¬ 
quency. 


D. Analytic continuation 


FIG. 1. Comparison of the error of the real part of the 
self-energy calculated with the conventional (red) and the im¬ 
proved (blue) method. Parameters are U = 20, fj, — 0.2317, 
T — 1. The improved method significantly reduces the error 
by making it linearly, rather than quadratically, dependent 
on frequency. 


C. Efficient calculation of self-energies 


The self-energies can be calculated by solving the 
single-impurity model ([^ with the CT-QMC method and 
then using Eq. (|^. However, this calculation is subject 
to stochastic uncertainties at high positive and negative 
frequencies for the following reason: In the first step the 
impurity problem is solved with the CT-QMC from which 
one obtains the Green function G(ia;„), which has the 
form G(iw„) = 5^ + 0{-^r^) for large |a;„|. The nu¬ 
merical error of this result is almost independent of a;„. 
Therefore relative error grows linearly with a;„. In the 
second step the Green function is inverted and substi¬ 
tuted into Eq. ([^. The linear term cancels out, and only 
terms of the order of (ia;„)®, (iw„)“^ and smaller remain. 
This results in a progressively larger relative error, i.e., 
the error increases quadratically with frequency (see 
Fig.0. 

To improve the accuracy of the calculation of the self¬ 
energies we adapt the method proposed by Snoek and 
Hofstetter,!^ and use the two-particle Green functions. 
The strategy goes as follows: Finding equations for the 
self-energies requires solving the equations of motion 
for the Green functions in the bosonic single-impurity 
model d^. This is similar to the method proposed for 
fermionsP^ The final expression for the self-energy reads 


(iw„) = 


[/F(icj„) + {k + fJ,) 1 ^nO 


C-\ 


( 11 ) 

where F(ia;„) is the Fourier transform of a matrix of two- 


The main goal of this paper is to calculate the mo¬ 
mentum resolved spectral function 7l(k, oj) and the mo¬ 
mentum integrated spectral function 7l(a;), respectively, 
of the Bose-Hubbard model Q. From the CT-QMG im¬ 
purity solver we obtain the Green functions and the self¬ 
energies in imaginary time or in Matsubara frequencies. 
Since these functions are analytic in the upper complex 
plane we can analytically continue them to the real axis. 
This involves an inversion of the Hilbert transform: 



It is well-known that this is a numerically ill-posed prob¬ 
lem, in particular for noisy data obtained with GT- 
QMG. To make the analytic continuation tractable we 
use the maximum entropy (MaxEnt) method.l^^M^ We 
found that the MaxEnt procedure sometimes fails to 
converg^^, especially in calculations of sharply peaked 
momentum resolved spectral functions A(k, w). In these 
cases we use the historic rather than the Bryan version 
of MaxEnt .1^ 


III. RESULTS 


In the following we present results for spectral func¬ 
tions of the Bose-Hubbard model on a simple cubic lattice 
for weak, intermediate, and strong interaction strengths. 
Eor illustration we indicate the sets of parameters at 
which the calculations for strong and intermediate cou¬ 
pling, respectively, were performed in the phase diagram 
(Fig. The diagram was calculated within the static 
Eisher mean-field theorj®! and the B-DMFTj ^°F0 | which 
was solved by GT-QMdsSM! and LGE.lial 
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FIG. 2. Phase diagram of the Bose-Hubbard model on the 
simple cubic lattice obtained with the static Fisher mean-field 
(MF) theorjl^and the B-DMFT, which is solved by CT-QMC 
and LCE. Only the first lobe corresponding to (n) ~ 1 is plot¬ 
ted. Circles, diamonds and triangles represent sets of parame¬ 
ters at which calculations for strong and intermediate interac¬ 
tions, respectively, were performed. Circles (orange): Strong¬ 
coupling regime; triangles (violet) and diamonds (blue): In¬ 
termediate interaction regime. 


A. Weak-coupling limit 


Our results in the weak interaction limit were obtained 
for T = 0.5, U = 0.25, and /r = —2.875. For these param¬ 
eters the average occupation per site is (n) = 0.4968 ± 
0.0005 and the average number of condensed bosons per 
site is (6)2 = 0.4542 ± 0.0004. Since U{n) « 0.1242 is 
small in comparison to the bandwidth 2zt = 6, the weak- 
coupling Bogoliubov approximatiorPi®^ can be expected 
to be applicable. The momentum resolved spectral func¬ 
tion A(k, w) obtained within the B-DMFT is presented 
in Fig. 1^. There are two bands: with positive energies 
for particle addition and with negative energies for parti¬ 
cle removal. Most of the spectral weight is concentrated 
in the npper band. The widths of the peaks for specihc 
k-points represent the mean lifetime of the quasiparti¬ 
cles. The weight of the lower band is orders of magnitude 
smaller, and therefore its exact shape and position is not 
determined reliably by MaxEnt. The dispersion relation 
e(k), obtained from A(k, w) according to the definition 

e(k) = max A(k, w), (15) 

UJ 

is shown in Fig. The dispersion relation obtained 
from the Bogoliubov approximation is found to be in very 
good agreement with the B-DMFT result, except near 
the F point, where the B-DMFT dispersion does not go to 
zero. This is attributed to the fact that the Hugenholtz- 
Pines theoreirP^is not fulfilled in the B-DMFT as already 
reported by Anders et al^ 


a) log(|A(k,(a)|) 
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FIG. 3. (a) Momentum resolved spectral functions at T = 0.5, 
U — 0.25, and fj, = —2.875 along the symmetry lines in the 
first Brillouin zone of a simple cubic lattice, (b) Results ob¬ 
tained by B-DMFT (red) and by the Bogoliubov approxima¬ 
tion (black) for the dispersion relations along the F-R line. 


B. Strong-coupling limit 

Calculations in the limit of strong interactions were 
performed for T = 0.5, and U = 50. In Fig. we 
present the spectral function A(k, w) for the Mott in¬ 
sulating phase at /r = 0.4C7 and the superfluid phase at 
/r = 0.9t7 (see Fig. [^. A striking difference between 
the spectra in the Mott insulating and superfluid phases, 
respectively, is the number of bands and their width at 
specific k-points. 

In the Mott insulating phase there are two bands, sep¬ 
arated by a gap whose width is approximately given by 
U — 3zt. This value becomes exact in the atomic limit 
(t —)■ 0). The factor izt corresponds to the sum of half¬ 
widths of the upper and lower bands. The shape of this 
dispersion relation is almost the same as that for a non¬ 
interacting band. In the Mott insulating phase with av¬ 
erage occupation (n) « 1 one particle is frozen on every 
lattice site such that an extra particle or hole can move 
almost freely through the system. Since the energy for 
creating an excitation is large, the density of excitations 
is very low, and therefore one can neglect their interac¬ 
tion. As a result quasiparticle and quasihole excitations 














5 


a) log(|A(k,(n)|);n=0.4U; Ml 
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-2.5 FIG. 5. Dispersion relation as calculated within the B- 

_3 DMFT and the static Fisher mean-field (MF) theorjl^, re¬ 

spectively. Top panel: Mott insulating phase at /i = 0.417 
and U = 50; bottom panel: Superfiuid phase at ^ = 0.917 
and U = 50. Only the three dominant bands are plotted. 


FIG. 4. Momentum resolved spectral function obtained from 
analytic continuation with MaxEnt of the GT-QMC data. 
Top panel: Mott insulating phase a.t — 0.417 and U = 50, 
bottom panel: Superfiuid phase at /r = 0.91/ and U = 50. 
Spectral functions are plotted along the symmetry lines in 
the first Brillouin zone for a simple cubic lattice. 

behave almost like noninteracting particles, i.e., have al¬ 
most noninteracting dispersion relations. The width of 
the peaks is therefore small, implying that the mean life¬ 
time is large. In fact, we checked that the width cor¬ 
responds to the uncertainty of the analytic continuation 
rather than to the mean lifetime. 

The bandwidth of particle excitations is wider than 
that of hole excitations by approximately 2zt. In the 
strong interaction limit the bandwidths of the hole and 
particle excitations in the Mott insulating phase with 
integer filling (n) are 2zt{n) and 2zt{{n) + 1), respec¬ 
tively. This is a quantum effect related to particle indis- 
tinguishability, which is simple to derive by starting with 
a state in which each site is occupied by n particles and 
treating the hopping term in the Hamiltonian Q as a 
small perturbation. 

In the superfiuid phase this is no longer valid. As the 
chemical potential increases, the energy for creating par¬ 
ticle excitations decreases. Therefore the interaction be¬ 
tween the quasiparticles needs to be taken into account. 
Indeed the spectrum of the superfiuid is significantly dif¬ 
ferent from that in the Mott insulating phase. In partic¬ 
ular, we find four rather than two bands. Only the po¬ 


sitions of the three lowest bands are determined reliably 
by MaxEnt. The four bands arise from two processes: 
(i) the splitting of a single band due to the interaction, 
similar to the Mott insulating phase; (ii) the creation of 
Bogoliubov quasiparticles due to the mixing of particle 
and hole excitations. We checked that, in contrast to the 
Mott insulating phase, the width of the peaks for small 
|w| is robust with respect to MaxEnt parameters as well 
as re-sampling and therefore represents the mean lifetime 
of quasiparticles and not the accuracy of analytic contin¬ 
uation. 

In Fig. we show the dispersion relations e(k) ob¬ 
tained within the B-DMFT scheme according to Eq. 
( |15[ ). We compare them with the dispersion relations ob¬ 
tained from the self-energies calculated within the static 
Fisher mean-field approximation.^ The results of both 
methods are in good agreement in the Mott insulating 
phase (Fig. top panel). 

In the superfiuid phase (Fig. bottom panel) this 
comparison is presented for the three dominant bands. 
The results are also in good agreement with the static 
Fisher mean-field results. In both approaches the high- 
energy (negative) band is similar to the band of a nonin¬ 
teracting hole. The remaining two dispersions are linear 
for small values of k = 0 in the Fisher mean-field theory 
and correspond to massless Bogoliubov quasiparticles. In 
the B-DMFT we also see the linear behaviour except 
for the vicinity of the T-point. The dispersion relation 
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FIG. 6. Spectral function A{uj) obtained with LCE (dotted 
line) and CT-QMC (full lines) at T = 1. Fade approximants 
were used for the data obtained with the LCE. Left panel: 
Calculations at ^ = 0.41/ for different values of U ; right panel: 
calculations at U = 15.5 for different values of fi. For U = 
12.5, fi = 0.4U and U = 15.5, ^ = 0.21/ the system is in the 
superfluid phase; for reference see Fig. 


around w = 0 does not go to zero, since the Hugenholtz- 
Pines theorem is not obeyed in the B-DMFT (see also 
Section III A). The slight deviation of the CT-QMC re¬ 
sult from that of the static Fisher mean-field theorj®! jn 
the middle of the band around the R point is attributed 
to the finite resolution of MaxEnt. 


C. Intermediate interaction and comparison with 
the strong-coupling solver 

The most interesting regime is that of intermediate 
interactions where both the Bogoliubov approximation 
and the static Fisher mean-field theorj®! are no longer 
valid. We computed the spectral functions in this regime 
within B-DMFT and compare our CT-QMC results with 
those obtained within the LCEp^. 

We used the same set of parameters as in Ref. [m 
i.e., each set corresponds to a point in a {^/U,t/U} pa¬ 
rameter space for T = 1. The selected points allow us 
to study the evolution of the spectral functions through 
the phase transition. Here we consider the phase transi¬ 
tions driven by the change of the interaction and by the 
change of the chemical potential (Fig left and right 
panel, respectively). 

In Figj^we present the B-DMFT results and compare 
with the LCE. There is good agreement regarding the 
widths and the positions of the bands in the superfluid 
phase. However, in the Mott insulating phase the LCE 
spectral functions develop a two-peak structure for both 


a) log(|A(k,co)l); U=18; Ml 
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FIG. 7. Momentum resolved spectral function A(k, ui) in the 
first Brillouin zone of a simple cubic lattice at ^ = 0.41/ and 
T = 1. Interaction strengths are the same as in Fig. left 
panel. 


positive nad negative energies. This is not supported 
by the CT-QMC results. The origin of this feature in 
the LCE may be an overfitting of the numerical analytic 
continuation, since in the LCE the error is small but un¬ 
known. 

We now focus on the momentum dependence of the 
spectral functions and the dispersion relations for the 
Mott insulating and superfluid phases (Figs. [^[^. The 
most striking feature is the difference of the width of the 
peaks for specific k-points in these two phases. Deep in 
the Mott insulating phase the width represents the un¬ 
certainty of the MaxEnt fit. This is not the case in the 
superfluid phase in which the width represents the mean 
lifetime of quasiparticles. Effectively, deep in the insu- 
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FIG. 9. (a): Momentum resolved spectral function at T = 0.5, 
U = 5, and fi — —2.625 along the symmetry lines in the first 
Brillouin zone of a simple cubic lattice, (b): Results obtained 
by the B-DMFT (red) and by the Bogoliubov approximation 
(black) for the dispersion relations along the F-R line. 


FIG. 8. Dispersion relation of correlated bosons for inter¬ 
mediate interaction strengths obtained within the B-DMFT. 
The momentum resolved spectral function 2l(k, uj) presented 
in Fig. [^was employed together with Eq. ( |15[ |. The results 
are compared to those of the static Fisher mean-field (MF) 
theorjP^ except for the value U = 14, where the B-DMFT 
finds the system to be in the Mott insulating phase, while it 
is in the superfluid phase according to the Fisher mean-field 
theory. 

lating phase we can describe the particles as almost free, 
whereas approaching the phase transition the quasipar¬ 
ticles obtain a mean lifetime. 

In the superfluid phase we observe high-energy excita¬ 
tions, which are shown in Fig. [^, resembling a structure 
found in the strong-coupling regime. Within the accu¬ 
racy of our method a more precise calculation of their 
position is not possible. 

The dispersion relations shown in Figj^ are compared 
with the strong-coupling dispersions obtained for [/ = 18 
and U = 12.5. In the Mott insulating phase the gap ob¬ 
tained by the static Fisher mean-field theorj®! jg much 


smaller than in the B-DMFT. Thereby the phase tran¬ 
sition is shifted to larger interactions, see Fig|^ In the 
superfluid phase we observe that, in contrast to the weak 
and strong-coupling limits, the dispersion relation ob¬ 
tained at intermediate coupling within the B-DMFT is 
not symmetric with respect to zero energy.l^ 

In Fig. we present the results obtained in the dilute 
gas regime. The parameters are T = 0.5, U = 5, and 
= —2.625. This corresponds to (n) = 0.12713±0.00001 
and (6) = 0.30301±0.00001. As in the previous cases (su¬ 
perfluid phase) we observe that the peaks in the spectral 
function are wide, which means that the mean lifetime 
of excitations is finite. The negative energy band has a 
weight which is orders of magnitude smaller, and there¬ 
fore its exact shape and position is again not determined 
reliably by our method of analytic continuation. The re¬ 
sults at the F-point are not included, since MaxEnt failed 
to converge for very small k. 

Although the particle density is low and the majority 
of the particles constitute the condensate, the interac¬ 
tion is not weak. Therefore one is not in the parameter 
range where the Bogoliubov approximation is applicable. 
Indeed, as seen in Fig. Ep. the dispersion relation is not 




































reproduced by this approximation. A similar discrepancy 
was reported in the experimeniP, where at high momenta 
the Bogoliubov approximation was found to overestimate 
the measured excitation energy. A further effect of the 
interaction is the appearance of an additional band at 
high energies, Fig. This high energy band appears 
to repel the one with low energies. As a result the low 
energy band is narrower. This explains the discrepancy 
between the B-DMFT and the Bogoliubov results, since 
the high energy band is absent in the latter. 


IV. SUMMARY 


where aj = (a| ai) and are vectors in Nambu nota¬ 
tion. Here ai (dj) annihilates (creates) a particle from the 
bath state I, b (&^) annihilates (creates) a particle on the 
impurity, V; is the Nambu matrix of couplings between 
the impurity and a bath state I, e; is the energy of the 
bath state I, and the vector represents the condensate 
field to which the impurity is coupled. The bath state 
energies ei and couplings V; are chosen such that 

A(ict;n) = ^ V/(ia;n<[73 - (A2) 

i 

where A(ia;„) is the dynamical mean field from Eq. ([^. 


We presented a computational method which allows 
one to calculate spectral functions for bosonic systems 
described with the Bose-Hubbard Hamiltonian. This ap¬ 
proach has the great advantage of being applicable for ar¬ 
bitrary values oiU/t, the ratio of the interaction strength 
and the hopping amplitude. It reproduces well the re¬ 
sults in the limiting cases of large and small values of 
the interaction. Results obtained by this method were 
shown for parameters where perturbative methods fail, 
thus providing new insights into the properties of corre¬ 
lated lattice bosons in previously inaccessible parameter 
ranges. The results show that while in the insulating 
phase particle and hole excitations tend to behave like 
free particles, in the superfluid phase the system is de¬ 
scribed by Bogoliubov quasiparticles with a finite mean 
lifetime. Finally, we discussed an alternative method for 
calculating the local self-energy in the CT-QMC solver 
which employs two-particle Green functions, and which 
produces accurate data at large Matsubara frequencies. 
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Appendix A: Calculation of bosonic self-energy from 
equations of motion 


We consider the Hamiltonian representatiorl^ of the 
impurity model ([^ 


H,mp= -tib^b+Wb^bb+ 


El 


+ eid\di^ — 


(AI) 


For this Hamiltonian we consider the disconnected 
Green functions 

G^oW = -(7^rb(r)bt(0)), (A3) 

Gfo(r)=-(TVaKr)bt(0)). (A4) 

Here the indices 0 and I refer to the impurity and bath 
states, respectively. The relation between Gqq (t) and the 
connected impurity Green function Goo(t) defined in Q 
is given by 

G^oM = Goo(T)-(b(r))(bt(0)). (A5) 


The time evolution of the operators is governed by the 
Hamiltonian Himp (for brevity we will drop the index 
henceforth), and we can calculate the derivatives of the 
Green functions with respect to imaginary time. For the 
local Green function we have 


a.ei!,(r) = -«(r)„3 - ^ j (*'(0) im) 'j 

= -d(r)(iJ3 -f/ xo3G(5o(t) - U(IT 3 F(t) 

-P3E^V^Gfo(r)-K03’®'^^ 

(A 6 ) 

where [• • • , • • • ] denotes the commutator, (IJ 3 is the diago¬ 
nal Pauli matrix, = (( 6 ^) (&)) is a Nambu vector, and 
the Nambu matrix F is defined as 


F(t) = 


- T. 


b^ (t)6(t)6(t)6'I’ (0) b^ (T)6(r)6(T)6(0) 
6'^ ('r)Sf (T)S(r)S'^ (0) ('r)Sf (t)&(t)S(0) 


a similar relation holds for G/qI 


drCfoir) = - 


[H,di]{T) 
[H,d]]{T) 

= -eidJsGfo - iiJ3Vj'G((o- 



(A7) 


(A 8 ) 


To handle the imaginary time derivative we perform a 


9 


Fourier transform to Matsubara frequencies and obtain 
ia;n(iT3G^o(i‘^n) = 1 - MGoo(iWn) + J2i '^iGfp(iw„) 


at 0 for non-zero Matsubara frequencies (for n ^ 0 
one has Gqq = Gqo = G). 


+ + Uf{iuin), 


(A9) 


0 : 


By inserting (All) into (A9) one finds in the case n = 


■ pdr ^ pdr ^ tAim = + + + 

ia;„(U3Gfo(ia;„) = e/Gfo(ia;„)-I-V]Ggo(ia;„). (AlU) ^ i 

_ ^ (A14) 

Using (AlO) it is easy to find an expression for Gfp in Starting from (|A2[) for n = 0 and replacing the discon- 


terms of Gqq: 

Cfo(i^«) = (iw„W 3 - e/1)"^ Vj'Goo(iw„). (All) 

We concentrate now on the case n 7 ^ 0, i.e., nonzero 
Matsubara frequencies. Inserting (All) into (A9) one 
obtains 

Goo(i^ra) = 

iWn<iT 3 + ^ Vi(iw„(U 3 - eil)"^vj' (1 -h C/F(iw„)). 

(A12) 


Using 


one can write 


Goo(iw„) = (ia;„<iT 3 -|-/rl-A(ia;„)) ^( 1 -|-C/F(ia;„)). (A13) 


Combining (A13) with the Dyson equation (^ 8 ^ we arrive 


nected Green function with the connected one leads to 

Goo(O) - = (Ail - A(0))-i(l -h + Uf{0)). 

(A15) 

In the last step the self-consistency condi tion ([^, i.e., 
= (k1 -I-A(0))$, is inserted into (A15). After some 
regrouping one obtains 

Goo( 0 )-$$^ = 

(Ail - A(0))-^[1 + {k + Ai)$$^ + C/F(0)] - 

(A16) 

After the above transformations one finally finds 

Goo(O) = (Atl + A(0))-i[l + (ic + Ai)$$^ + UF(0)], (A17) 

which, in combination with the Dyson equation (|^, re¬ 
duces to 0 for n = 0 . 
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